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^ ■ Abstract 

fj . We present simple equations for a canonical-basis formulation of the time-dependent Hartree- 

q^ I Fock-Bogoliubov (TDHFB) theory. The equations are obtained from the TDHFB theory with 



p 

o 
o 



an approximation that the pair potential is assumed to be diagonal in the canonical basis. The 
canonical-basis formulation significantly reduces the computational cost. We apply the method to 



> 

in 

00 . 

r~^ , linear-response calculations for even-even light nuclei and demonstrate its capability and accuracy 



by comparing our results with recent calculations of the quasi-particle random-phase approximation 
with Skyrme functionals. We show systematic studies of El strength distributions for Ne and Mg 
isotopes. The evolution of the low-lying pygmy strength seems to be determined by the interplay of 
several factors, including the neutron excess, separation energy, neutron shell effects, deformation. 



>< 
C^ • and pairing. 



I. INTRODUCTION 

The time-dependent Hartree-Fock (TDHF) theory was extensively apphed to studies of 
nuclear collective phenomena as a microscopic approach to nuclear dynamics[l|. Recently, 
it has been revisited with modern energy density functionals and more accurate descrip- 
tion of nuclear properties has been achieved 2H,7l] . The TDHF theory uses only occupied 
orbitals, number of which is equal to the number of particles (A^), to describe a variety 
of nuclear dynamics such as heavy-ion scattering, fusion/fission phenomena, and linear re- 
sponse functions. However, it neglects the residual interactions in particle-particle and 
hole-hole channels, which becomes problematic especially for open-shell heavy nuclei. An 
alternative approach including the pairing correlations is the time-dependent Hartree-Fock- 
Bogoliubov (TDHFB) theory 8|] . The TDHFB equation is formulated in a similar manner to 
the TDHF, however it uses the quasi-particle orbitals instead of the occupied orbitals. Since 
the number of the quasi-particle orbitals is, in principle, infinite, the accurate calculation of 
TDHFB is presently impractical. Only recently, a few attempts of the TDHFB calculation 
have been performed, but either with a small model spacejot] or with restriction to spherical 



symmetry 



lOJ. 



A much simpler approach was proposed by Blocki and Flocard in Ref. [ll(]. They gave 
equations of motion for time-dependent canonical states \4>k{t)) {k = 1,- ■ ■ , M) and those 
for the time-dependent BCS factors {uk{t),Vk{t)). Since the number of canonical basis is 
larger than the particle number but not significantly different (M ~ A^), the necessary 



computational task is roughly the same as that of TDHF. The similar met 



lods have been 



- |l4J . However, 



applied to studies of heavy- ion reactions with use of simple functionals [12- 

it has never been tested with realistic modern functionals so far, and we do not know how 

reliable this approximated scheme is. In addition, although the equations of motion were 



provided for a very schematic pairing functional in Ref. ll| , its theoretical foundation seems 
rather obscure to us. 

In this paper, we derive the equations of motion for general functionals and clarify the 
approximations/assumptions we need to make. We call those equations "Canonical-basis 
TDHFB" (Cb-TDHFB) equations. We apply the method to the linear-response calculations 
using the full Skyrme functionals. The results will be compared with recent calculations 
of the quasi-particle random-phase approximation (QRPA), then demonstrate its feasibility 



and accuracy. 

The paper is organized as follows. In Sec. HIl we present the basic equations of the 
present method and their derivation. Especially, we would like to clarify what kind of 
assumption/approximation is necessary to justify the Cb-TDHFB equations. In Sec. IIIH we 
show properties of the Cb-TDHFB, including gauge invariance, conservation laws, and the 
small- amplitude limit. In numerical calculations in this paper, we adopt a schematic choice 
for the pairing energy functional, similar to Ref . [ll| , which violates the gauge invariance. In 
Sec. IIVI we show effect of the violation of the gauge invariance can be minimized by a special 
choice of the gauge condition. In Sec. El details of our numerical installation is given. Then, 
in Sec. [VI\ we present numerical results of the real-time calculations of the linear response 
and compare them with recent QRPA/RPA calculations. Finally, the conclusion is given in 
Sec.lVin 

II. DERIVATION OF BASIC EQUATIONS 

In this section, we derive the basic equations of Cb-TDHFB method. Using the 
time-de pen dent variational principle, the similar equations were derived by Blocki and 



FlocardpjJ. However, it was not clear that what kind of approximation was introduced 
and how they are different from the full TDHFB. We present here a sufficient condition to 
reduce the TDHFB equations to those in a simple canonical form. 

We start from the density-matrix equation of the TDHFB and find equations for the 
canonical-basis states and their occupation- and pair-probability factors. In order to clarify 
our heuristic strategy, let us start from a simpler case without the pairing correlation. 

A. TDHF equation 



The TDHF equation in the density- matrix formalism is written as 15 1 



^^pit) = [hit),pit)], (1) 

where p{t) and h{t) are the one-body density operator and the single-particle (Hartree- 
Fock) Hamiltonian, respectively. We now express the one-body density using the time- 
dependent canonical single-particle basis, {\4>k(t))}, which are assumed to be orthonormal 



i{Mt)\Mt)) = Ski). 



N 



p{t) = J2\Mt)){Mt)l (2) 

fc=i 

where N is the total particle number. Substituting this into Eq. ([1]), we have 

N N 

J]{^|0,(t))(0,(t)|+z|0,(t))(0fc(t)|} = J]{Mt)|0fc(t))(0,(t)|-|0fc(t))(0,(t)|/i(t)}. (3) 
fc=i fc=i 

the inner product with \(pk(t)) leads to 

P (i^ - h{t)] \Mt)) = k = l,---,N, (4) 

with P = 1 — J2k=i\'^k(t)){(pk(t)\. Here, we used the conservation of the orthonormal 
property for the canonical states, d / dt{(l)k{t)\(l)i{t)) = 0. This leads to the most general 
canonical-basis TDHF (Cb-TDHF) equations 

TV 
, \rh,Jt\\ = h(tMrh,Jt\\ -' 

dt 



1=1 



where the matrix riik{t) is arbitrary but should be hermitian to conserve the orthonormal 
property. It is easy to see that the time evolution of the density does not depend on the choice 
of r]ik. This is related to the gauge invariance with respect to the unitary transformations 
among \4>k{t)) {k = 1, ■ ■ ■ ,N). The most common choice is rjik = 0, which leads to the 
TDHF equation shown in most textbooks. 

B. Cb-TDHFB equations 

We now derive Cb-TDHFB equations starting from the generalized density-matrix for- 
malism. The TDHFB equation can be written in terms of the generalized density matrix|8| 



as 



where 



i^R=[H,R], (6) 



-K* l-p* \-A* -h* 



This is equivalent to the following equations for one-body density matrix, p{t), and the 
pairing-tensor matrix, K,{t). 

4-^^ = [^W'PW] + <tWit) - ^it)^*it), (8) 

z^n{t) = h{t)K{t) + K{t)h*{t) + A(t)(l - p*{t)) - p{t)A{t). (9) 

At each instant of time, we may diagonalize the density operator p in the orthonormal 
canonical basis, {(j)k{t), 4>k{t)} with the occupation probabilities pk- For the canonical states, 
we use the alphabetic indexes such as k for half of the total space indicated by A; > 0. For 
each state with k > 0, there exists a "paired" state A; < which is orthogonal to all the 
states with k > 0. The set of states {(j)k,4'k} generate the whole single-particle space 28']. 
We use the Greek letters /i, z/, ■ ■ ■ for indexes of an adopted representation (complete set) for 
the single-particle states. The creation operator of particles at the state \(pk(t)) is expressed 
as c^(t) = J2n{f^\^ki't))cjj^, and the TDHFB state is expressed in the canonical (BCS) form 
as 

i^w) = n {^fcW + MMiHit)} io). (10) 

fc>0 

For later purposes, it is convenient to introduce the following notations for two-particle 
states: 

{pu\Mt)Mt)) = {p\Mt)){i^\Mt)), (11) 

{{pu\Mt)Mt))) = {pi^\Mt)Mt)) - {pi^\4>-k{t)Mt)), (12) 

and for projection operator on a canonical pair of states {k, k), 

T^kit) ^ |0fc(t))(0.(t)| + |0fc(t))(0s(t)|. (13) 

Then, it is easy to show the following properties (/c, Z > 0): 

Y,{H<t>k{t)(t>-k{t)){Ut)Ut)\^^) = Ski, (14) 

J2{{HMt)Mm{{MtMt)\piy)) = 25ku (15) 

5^((/ia|0fc(t)0fc(t)))((0Kt)0r(t)k^)) = 5ki{pHt)k\y), (16) 

5^((/i(T|0fc(t)0fc(t)))(z/|7r^|a) = 6ki{{py\<pk{t)<p-k{t))). (17) 



Using these notations, the density and the pairing-tensor matrixes are given by 

P^i) = $^Pfe(t)(/i|7rfc(t)|z/), (18) 

fc>0 

^^^At) = j2Kk{t){{fiu\Mt)4>-km, (19) 

fc>0 

where pk{t) and Kk(t) are occupation and pair probabihties, respectively. In terms of the 



BCS factors of ('U,f)[l5|, they are given as pk(t) = \vk{t)\^ and Kk(t) = ul{t)vk{t). It should 
be noted that the canonical pair of states, \(f)k{t)) and |0^(t)), are assumed to be orthonormal 
but not necessarily related with each other by the time reversal, |0j;) 7^ T|0fc). 

Thanks to the orthonormal property, we can invert Eqs. flTSj) and flT^ for pk and Kk, 

Pk{t) = J2^Mt)\p)p,.it){u\Mt)) = J2^Mt)\p)p,.it){i^\Mt)), (20) 

With help of Eq. (ITSll . the derivative of Pk(t) with respect to time t leads to 

fJ,V 



J2(Mt)\p)^^{'^\Mt)) 



IJ,U 



= J2 {^^k{t)A;,{t){iyp\Mt)Mt)) + 4{t)A,At){Mt)Mt)\P'y)} ■ (22) 

HI/ 

We used the assumption of norm conservation for the second equation, and used the TDHFB 
equation (|H]) in the last equation. Since the pair potential A^,^(t) is anti-symmetric, it is 
written in a simple form as 

^j/k{t) = '^k{t)Al{t) - <(t)A,(t), (23) 

Afe(t) = -5^A^,(t)(0fc(t)0s(t)|^z/) = -i^ A^,(t)((0fc(t)0s(t)|/iz/)). (24) 

fit/ jJiU 

In case that the pair potential is computed from a two-body interaction v as A^y{t) = 
y^ a v^y^api^ajsif)) the gap parameters, Afc(t), are identical to those of the BCS approximation 



15|- 



^k{t) = - ^ i^i{t){vkk,u - Vkkfl) = -Y^ i^i{t)Vkk,u- (25) 



/>0 l>0 



It should be noted here that the two-body matrix elements V/^j^ii (and the anti-symmetric 
■y^j^;;-) are time-dependent because the canonical basis, {k,k) and {1,1), are time-dependent. 
In the same way, we evaluate the time derivative of K,k{t) as 

^j^^kit) = J2(Mt)Mt)\H^^ + ^^k{t) U^lMt)) + i^l^dt))) ■ (26) 

Then, using the TDHFB equation (^, we obtain 

^J^f^kit) = Mt) + r/fc(t)) Kkit) + Afc(t) {2pk{t) - 1) , (27) 

where r,k{t) = {Mt)\Ht)\Mt)) + ^{^\Mt))- 

The time-dependent equations for pk{t) and K,k{t) are now given in rather simple forms as 
Eqs. ( 123|) and ( 1271) . So far, their derivation is solely based on the TDHFB equations, utilizing 
the fact that p^u{t) and K,^u{t) can be expressed by the orthonormal canonical basis, |0fc(t)) 
and |0fc(t)), and their occupation and pair probabilities, Pk{t) and K,k{t)- However, in general, 
the time evolution of the canonical basis is not given by a simple equation. Therefore, we 
now introduce an assumption (approximation) that the pair potential is written as 

A^,(t) = -Y,Mt){{HMt)Mt)))- (28) 

k>0 

This satisfies Eq. ( 1M1) . but in general, Eq. ( l24l) can not be inverted because the two-particle 
states |0fc0fc) do not span the whole space. In other words, we only take into account the pair 
potential of the "diagonal" parts in the canonical basis, A^;- = —AkSki- In the stationary 
limit (|0fc) = T|0fc)), this is equivalent to the ordinary BCS approximation (see Sec. IIIICI) . 
With the approximation of Eq. fl28l) . it is easy to see that the TDHFB equations, ([8]) 
and (IHl), are consistent with the following equations: 

^Q^iMt)) = m) - vk{t))\Mt)). ^^\<t>-kit)) = (hit) - v-kimMt))- (29) 

In summary, the Cb-TDHFB equations consists of Eqs. (12^ . (125]) . and fl7r|) . To derive 
these equations from the TDHFB equations, we have assumed the diagonal property of the 
pair potential, Eq. ( I28ll . 



III. PROPERTIES OF THE CB-TDHFB EQUATIONS 

A. Gauge invar iance 

The rjk{t) and ?7^(t), in Eqs. (1271) and fl2^ . must be real to conserve the orthonormal 
property, however, they are arbitrary. This is related to the phase degrees of freedom of the 
canonical states. The Cb-TDHFB equations, (!23|) . ( 127|) and (129|) . are invariant with respect 
to the following gauge transformations with arbitrary real functions, 9k{t) and 9^{t). 

10,) ^ e^^'^-«|0,) and |0s) ^ e^'^W|0s) (30) 

^, ^ e-^(^^W+'^^W)/€fc and A, ^ e-*(^^W+^sW)Afc (31) 



simultaneously with 



Vkit) ^ rik{t) + — and Vki^) ^ rj-^it) + —. 



The phase relations of Eq. ( 1311) are obtained from Eqs. ( 1211) and ( 1M1) . 



B. Conservation laws 

i. Orthonormality of canonical states 

Apparently, Eq. ( 129|) conserves the orthonormal property of canonical states, as far as rjk 
are real. 

^^{Mt)\Ht)) = (MmiHt) - vi{t)) - {h^{t)-Vk{t))}\Mt)) = 0. (32) 

Here, we assume (0fc(t)|0/(t)) = 6ki at time t. 

2. Average particle number 

The average particle number also conserves because 

^|iV(t) = 2^J^ Yl P>^^^) = 2 J2^'^kit)Alit) - 4(t) A,.(t)) = 0, (33) 

fc>0 A:>0 

where we used the expression of the pairing energy, Eq. (|59|) . for the last equation. 



8 



3. Average total energy 

Time variation of the energy functional E[p, k] can be divided into two: dE/dt = 
dE/dt\p + dE/dt\t^. The variation of energy associated with the normal- density fluctua- 
tion is 



. dE 
^ Itt 



where efc(t) = {(j)k{t)\h{t)\(f)k{t)) . This equation has an intuitive physical interpretation. The 
energy carried by a canonical state |0fc) is efc(t) x p^. If the occupation probability is fixed 
during the time evolution, the right-hand side of Eq. ( !34l) vanishes. This corresponds to 
the cases such as the TDHF and its extension with fixed BCS occupation probabilities. In 
the TDHFB, the energy variation in Eq. ( l34l) transfers from/to the pairing energy. In fact, 
time variation due to the pairing tensors-produce fluctuation produces 

dE 



I —r- 

dt 



^EE%^ + |J^) -E«^--^«(M<)+..W). (36) 



I 



where Eq. (|28|) is used. Because of Eq. (123|) . two contributions of Eqs. t^^ and ( 135|) always 
cancel and the total energy is conserved. This is natural because the Cb-TDHFB equations 
satisfy the TDHFB equations, (|8]) and ([9]), for which the conservation of the total energy in 
TDHFB is well known|8|. 



C. Stationary solution 

When we assume that all the canonical states are eigenstates of the time-independent 
single-particle Hamiltonian Hq. 

\Mt)) = \4)e'''^'\ \Mt)) = me''-^^'\ (36) 

ho\<Pl) = el\4), ho\<Pl) = el\<Pl), (37) 

where \(f)%) = T|{/)fc) have the same eigenvalues e° as \4>k)- Here, dOk/dt = i{d(j)k/dt\(j)k) 
and d9j,/dt = i{d(f)j,/dt\(l)j,) are arbitrary real functions of t. K,k{t) and Afc(t) should have 
a common time-dependent phase associated with the chemical potential A as e"^'"^*. In 
addition to this, according to their definition, Eqs. (I2T1) and ( l24l) . they have the following 



additional phases connected with the phases of the canonical states. 

Kkit) = 4exp{-i{2Xt + Okit) + e-^{t))}, (38) 

A,(t) = A° exp{-z(2At + O^it) + 0^(t))}, (39) 

The stationary case of Eq. (12^ . dp^/dt = 0, indicates that k^ and A° have the same 
arguments to make fi;fc(^)A^(t) real. If all the pairing matrix elements are real, we can 
choose both n^ and A^ are real. Then, Eq. (l27|l reduces to 

2{el-X)KUAU2pl -!) = (}. (40) 

This is consistent with the ordinary BCS result. 

pI = -\l- '° ~ ^ — 1 , (41) 



1 AO 



D. Small amplitude limit and the Nambu-Goldstone modes 



(42) 



It is known that the small-amplitude approximation for the TDHFB around the HFB 
ground state is identical to the QRPA. In the QRPA, when the ground state (HFB state) 
breaks continuous symmetry of the Hamiltonian, the Nambu-Goldstone modes appear as the 
zero-energy modes. In this section, we show that this is also true for the small-amplitude 
limit of the Cb-TDHFB. 

The ground state is given by |0°), |0^), k^, and p° which satisfy Eqs. ( l37|l and (HO|l . 
Extracting trivial phase factors, ikif) = /q {^fc(i') ~ ^2}'^^'; we express the time-dependent 
quantities as follows: 

\Ut)) = m + \Hk{t))] e^^'^W, \<p-,{t)) = {\4) + \6Mt))} e'^-^^'\ (43) 

n,{t) = {4 + 5«:fc(t)} e-^{«^W+«^W+2At}^ ^^(^) ^ |^o ^ 5Afc(t)} e-^{«^W+«^W+2"*>(,44) 

Pkit) = pl + 6pkit), h{t) = ho + 6h{t), (45) 

Substituting these into Eqs. (|29|1 . fl23l) . and fl27j) . they lead to the following equations in the 



10 



linear order with respect to the fluctuation. 

^§-^\SMt)) = {ho - 4)\SMt)) + mMl), (k ^ k), (46) 

tpPkit) = Al*6K,{t) + 4SAl{t) - c.c, (47) 

t^^6^k{t) = 2(e° - X)6Kk{t) + {2pl - l)5Afe(t) + 2A°5pfc(t). (48) 

When these fluctuating parts have specific oscillating frequency cj, they correspond to the 
normal modes. The zero-energy modes correspond to stationary normal-mode solutions with 
a; = 0. 



1. Translation and rotation 

When the HFB ground state spontaneously violates the translational (rotational) sym- 
metry, there are generators, P (J), which transform the ground state into a new state but 
keep the energy invariant. Here, let us denote one of those hermitian generators, S. The 
transformation with respect to the generator 5* with real parameter a leads to 

10^ ^ \<fM)) = e'^'\ct>l) (k^k), (49) 

ho -^ ho{a) = e^"^/ioe-^"^, (50) 

with pk{a) = Pfc, i^k{c() = /tfc, efc(a) = e°, and Ak{a) = A°. These transformed quantities 
should also satisfy Eq. (l37|l . 

{ho{a) - elMlia)) = 0. (51) 

In the linear order with respect to the parameter a, we have 

za{h - el)S\4) + za[S, hoM) = 0. (52) 

Equation fl52l) means that |50f ) = iaS\(f)1) and 6hs = ia[S, Hq] correspond to a normal- mode 
solution with w = for Eq. (jl6]). 5pf = 0, 6n^ = 0, and (5Af = also satisfy Eqs. (j47]) and 
( HHj) . Therefore, the Nambu-Goldstone modes related to the spontaneous breaking of the 
translational and rotational symmetries become zero-energy modes in the small-amplitude 
Cb-TDHFB equations. 



11 



2. Pairing rotation 

When the ground state is in the superfluid phase, we have k^ 7^ at least for a certain k. 
The ground state can be transformed into a new state by operation of e*^^ where N is the 
number operator. This transformation changes the phase of Kk and A^ but keeps the other 
quantities invariant. 

5k^ = e^''4 -4^ 2i94, SA^ = e^^'Al - A° ^ 2t9Al (53) 

6p^ = 5h^ = 0, |<50f ) = |50f ) = 0. (54) 

Using Eq. f HU]) . it is easy to see that these quantities correspond to an a; = solution of 
Eqs. (H^ . ( H7|l . and (HSj) . Thus, the pairing rotational modes appear as the zero energy 
modes as well. 

3. Particle-particle (hole-hole) RPA 

The Cb-TDHFB equation fH^ seems to be independent from the rest of equations, f l47p 
and ( 1481) . at first sight. However, this is not true in general, because 5Ak{t) depend on 
\^'Pk{t)) and 5h{t) depends on 6pk{t). In contrast, when the ground state is in the normal 
phase {k^ = A° = 0), SA^it) becomes independent from \6(f)i(t)), and we have Spk{t) = 0. 
This means that the particle-hole (p-h) channel is exactly decoupled from the particle- 
particle (p-p) and hole-hole (h-h) channels. It is well-known that TDHF e quat ions in the 
small-amplitude limit, (HH]), reduce to the RPA equation in the ph-channeljsl. llSl Il6|. Thus, 
we here discuss properties of the p-p and h-h channels. 

The p-p and h-h dynamics are described by the following equations. 

t^6Kk{t) = 2el6Kk{t) ± 6Ak{t), (55) 

where the sign + (— ) is for hole (particle) orbitals, and we omit the chemical potential 
A. For the p-p channel {u = En+2 — E^), a normal mode with frequency u is described 
by Sup = XpC"*"^* for particle orbitals {\p\ > N/2) and Suh = — y/iC^"^* for hole orbitals 
{\h\ < N/2). For the h-h channel {uj = En-2 — En), it is described by 6Kh = XhC'"^^^ for 
hole orbitals {\h\ < N/2) and 6Kp = — Y^e^*"^* for particle orbitals (|A;| > N/2). Equation 



12 



(155|) can be rewritten in a matrix form as 

-Vhhp'p' -"^^IShh' + Vhhh'h' J \^h'j \o -I J yZh^ 

where Zp = Xp {Zp = Yp) and Zh = Y^ {Z^ = X^) for the p-p (h-h) channel. This is 



equivalent to the p-p and h-h RPA in the BCS approximation 15 1. 



IV. CB-TDHFB EQUATIONS WITH A SIMPLE PAIRING ENERGY FUNC- 
TIONAL AND GAUGE CONDITION 

A. Pairing energy functional 

Normally, the pairing energy functional is bi-linear with respect to K^jy and k* . For 
instance, when it is calculated from the two-body interaction, it is given by 

^«(^) = 5Z ^t^'^.P^'^lu{t)>^pcT{t)- (57) 

Thus, the pairing energy can be also written as 



= -J2 ^k{t)Am = -J2 4(i)Afc(t). (59) 

fc>0 A:>0 

For numerical calculations in the present study, we adopt a schematic pairing functional 
in a form of 

Egit) = -Yl GMKl{t)Ki{t),= -Y,4{t)Mt), Afc(t) = J^GkMt), (60) 

k,l>0 k>0 l>0 

where Gm is a hermitian matrix. This pairing functional produces a pair potential which 
is diagonal in the canonical basis. This is consistent with the approximation of Eq. fl28l) . 
However, the functional violates the gauge invariance (l3T|l . because 

l>0 k>0 

The violation comes from the fact that the Afc(t) in this schematic definition no longer hold 
the correct phase relation to canonical states {k, k), according to the definition of Eq. (I24l) . 

13 



Therefore, we require the gauge condition of (-^|0fc) = {-^\4>k) = 0, so as to minimize the 
phase change of canonical states. This means that we choose the gauge parameters rik{t) as 

Vk{t) = efc(t) = {Mt)\h{t)\Mt)), Vkit) = e-,{t) = (</)g(t)|/i(t)|</>s(t)). (62) 

B. Properties of Cb-TDHFB equations with Eg 

The Cb-TDHFB equations with the simple pairing functional ( 160|) keep the following 
desired properties, if we adopt the special gauge condition f l62|) . The details are presented 
in Appendix. 

1. Conservation law 

(a) Conservation of orthonormal property of the canonical states 

(b) Conservation of average particle number 

(c) Conservation of average total energy 

2. The stationary solution corresponds to the HF+BCS solution. 

3. Small-amplitude limit 

(a) The Nambu-Goldstone modes are zero-energy normal-mode solutions. 

(b) If the ground state is in the normal phase, the equations are identical to the 
particle-hole, particle-particle, and hole-hole RPA with the BCS approximation. 

Among these properties, 1(a) and 1(b) do not depend on the choice of the gauge, however, 
the other properties are guaranteed only with the special choice of gauge (1621) . 

V. DETAILS OF NUMERICAL CALCULATIONS 

A. Treatment of the pairing energy functional 

In numerical calculations, we start from the HF+BCS calculation for the ground state. 
The pairing energy is calculated for the constant monopole pairing interaction with a smooth 
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truncation for the model space. We follow the prescription given by Tajinia et al IT], which 
is equivalent to the following choice of Gki of Eq. fl60|) . 



Gu = 9f{el)f{e',), 



(63) 



with a constant real parameter g. The cut-off function /(e), which depends on the ground- 
state single-particle energies, is in the following form 



f{e) 



1 + exp 



e — tr 



0.5 MeV 



-1/2 



9{ec-e), 



(64) 



with the cut-off energies 



A + 5.0 MeV, Cc 



2.3 MeV, 



(65) 



where A is the average of the highest occupied level and the lowest unoccupied level in the 
HF state. Here, the cut-off parameter Cc is necessary to prevent occupation of spatially 
unlocalized single-particle states, known as the problem of unphysical gas near the drip line. 
For neutrons, if Cc becomes positive, we replace it by zero. 

To determine the pairing strength constant g for each nuclei, we again follow the pre- 
scription of Ref . 17| which is practically identical to the one in Ref . 18| . For light nuclei 
{A < 50), we replace g with 0.6 MeV when the calculated value exceeds 0.6 MeV. The 
pairing force strengths Gki are calculated for the ground state and kept constant during the 
time evolution. We define the state-independent pairing gap as follows: 



A(t)=^5^«:.(t)/(e°). 



(66) 



fc>0 



The gap parameter for each canonical pair of states, k and fc, can be written as Ak(t) 



B. Energy density functional and coordinate-space representation 

In the present calculations, we adopt a Skyrme energy functional, i?sky[p], with the pa- 
rameter set of SkM* 191] . The functional contains both time-even and time-odd densities 
same as Ref. 20|. The pairing energy functional is added to this, to give the total energy 
functional, E[p, k] = Es-ky[p] + Egln]. 
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We use the Cartesian coordinate-space representation for the canonical states, 
0fc(r, cr;t) = {f, a\(j)k{t)) with a = ±1/2. The three-dimensional (3D) coordinate space 
is discretized in square mesh of Ax = Ay = Az = 0.8 fm in a sphere with radius of 12 fm. 
Thus, each canonical state is represented by 0fc(z, j. A;, a; t) with three discrete indexes for 
the 3D space. 

C. Calculation for the ground state 

First, we need to obtain a static solution to construct an initial state for the time- 
dependent calculation. The numerical procedure is as follows: 

1. Solve Eq. f p7|) for occupied canonical states {\k\ < N/2) with Ok = 1, and construct 



the HF Hamiltonian hQ[p\, using the imaginary-time method 21 1. 



2. Calculate unoccupied canonical states (^°(r, a) (|fc| > N/2) below the energy cut-off 
Cc, using the imaginary-time method with Hq. 

3. Solve the BCS equations jl 51] to obtain pk and k^. 

4. Update /io[p] with new pk, then solve Eq. ( 137|) again with the imaginary-time method, 
to calculate canonical states with e^ < Cc. 

5. Back to 3. and repeat until convergence. 



To solve Eq. ( 1371) . the imaginary-time-evolution operator for a small time interval At is 
repeatedly operated on each single-particle wave function. After each evolution, the single- 
particle wave functions are orthonormalized with the Gram-Schmidt method from low- to 
high-energy states. We add the constraints for the center of mass, jfp{f) = 0, and the 
principal axis, J rirjp{r) = {i ^ j) for deformed nuclei. 

D. Real-time calculation for strength functions 

The canonical states in the ground state define the initial state for time evolution. In 
order to study the linear response, we use an weak instantaneous external field of Kxt(^; t) = 
—rjF{f)5{t) to start the time evolution. Here, F{f) is a one-body field, such as El operator 
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with recoil charges, 

{Ne/A)ri for protons 
F^{^ = r , (67) 

— {Ze/A)ri for neutrons 

where i = {x,y,z). We also study the isoscalar quadrupole response with 

F{r) = , \ ^ {r'Y2K{r) + r^Y^^r)} , 7^ = and 2. (68) 

'2(1 + 6x0) 



Then, at time t = 0+, the canonical states are given by 

Mr, a; t = 0+) = e'^^^^^r, a), 0s(f, a; t = 0+) = e'^^^^'Ulir, cr), (69) 

and the BCS factors by 

p,(t = 0+)=p°, Kfc(t = 0+) = K°. (70) 

The parameter rj controls the strength of the external field. In this paper, since we calculate 
the linear response, it should be small enough to validate the linearity. 

To solve the Cb-TDHFB equations in real time, we use the simple Euler's algorithm. 

i</)fc(t + 2dt) = i(j)kit) + {h{t + dt) - ek(t + dt)}(t>k{t + dt) x 2dt, (71) 

ipkit + 2dt) = ipkit) + {Kkit + dt)/\l{t + dt) - c.c.} x 2dt, (72) 
iKkit + 2dt) = iKk{t) + [Kk{t + dt){tk{t + dt) + ei{t + dt) - 2 A} 

+Afe(t + dt){2pk{t + dt)-l}]x 2dt. (73) 



Here, we insert the chemical potential in Eq. fl27j) which cancels a global time-dependent 
phase at the ground state, e"^*'^*, for k^ and A^. To construct the states at the first step 
of t = dt, we use the fourth-order Taylor expansion of the time-evolution operator for the 
canonical states [2[| and use the Euler's method for pk{dt) and Kk{dt). The time step dt is 
0.0005 MeV-^ The time evolution is calculated up to T = 10 MeV'^ 

The strength function with respect to the operator F is calculated with the following 

H 

formula [2j. 

SiE;F)^y2\{^^\F\%)\-'5{E-K) = -—lmfiE), K>0, (74) 

^-^ 7171 

n ' 

where E^ = En — Eq and f{E) is defined by 

/•oo /• 

fiE) = J t/te(^^'r/2)*yF(f){p(f,t)-p(f,0)}rff, (75) 
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where we have introduced a smoothing parameter F which is set to 1 MeV throughout the 
calculations in Sec. I VII The formula can be obtained from the time-dependent perturbation 
theory in the first order with respect to r] [2|. Note that the strength function S{E;F) is 
independent from magnitude of the parameter rj as far as the linear approximation is valid. 
In the present study, we adopt the value oirj = 10~^ fm~^ for the El operator, and t] = 10~^ 
fm~^ for the quadrupole operator. 

VI. NUMERICAL RESULTS OF LINEAR RESPONSE CALCULATION 

In this paper, we apply the Cb-TDHFB method to calculation of the strength functions 
for Ne and Mg isotopes. First, we show, in Table HI calculated ground-state properties; 
deformations, chemical potentials, and gap energies defined by Eq. fl66l) . These nuclei show 
a variety of shapes (spherical, prolate, oblate, and triaxial), with and without superfluidity. 
For nuclei in the superfiuid phase with A 7^ 0, the numbers of canonical orbitals, Mr, 
included in the calculation (e^ < Cc) are as follows: For protons, Mp = 16 for 24,26, 28]\jg ^^^^ 
for 26,28,30,32^g^ g^^^ ^^ ^ 20 for 30Ne. For neutrons, M„ = 20 orbitals for ^^Ne, M„ = 24 for 
^^Ne, Mn = 28 for 30.34,36 ]y[g^ j^^ ^ 30 fQj, 38,40]y[g^ These numbers are, of course, larger than 
the proton and neutron numbers, but not significantly different. In the case with A = 0, 
of course, we only calculate the occupied orbitals (Mp = Z and M„ = A^). Therefore, the 
numerical task of the Cb-TDHFB is in the same order as that of the TDHF. Note that, in 
the real-time calculation, the numerical cost is proportional to Mn + Mp. 

A. Isoscalar quadrupole excitations: Comparison with QRPA calculations 

We expect that the strength functions calculated in the present real-time approaches 
reproduce those in the QRPA. This is strictly true if we solve the full TDHFB equations, 
however, since we solve the Cb-TDHFB equations with the schematic pairing functional of 
Eq. ( 160|) with (!63|) . let us first show the comparison between our result and the HFB-f-QRPA 
calculations. The QRPA calculations have been done with a computer program for axially 



deformed nuclei developed in Ref. [22|, which diagonalizes the QRPA matrix of large di- 
mensions in the quasi-particle basis. This is based on the HFB ground state calculated 
in the two-dimensional coordinate-space representation with the Skyrme functional SkM* 
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TABLE I: Calculated ground-state properties of Ne and Mg isotopes; quadrupole deformation 
parameters (/3,7), pairing gaps ([66]) for neutrons and protons (A„,Ap), chemical potentials for 
neutrons and protons (A„,Ap). In the case of normal phase (A = 0), we define the chemical 



potential as the single-particle energy of the highest occupied orbital, An = e^ and Ap 
pairing gaps and chemical potentials are given in units of MeV. 



€%. The 





/? 


7 


A„ Ap 


— An —Ap 


20Ne 


0.37 


0° 


0.0 0.0 


13.07 9.19 


22Ne 


0.37 


0° 


0.0 0.0 


11.03 12.38 


24Ne 


0.17 60° 


0.0 0.74 10.57 13.04 


26Ne 


0.0 


- 


0.0 1.00 


7.17 14.92 


28Ne 


0.0 


- 


0.79 1.01 


3.22 17.05 


30Ne 


0.0 


- 


0.0 1.01 


3.79 19.09 


32Ne 


0.36 


0° 


0.95 0.0 


2.16 23.61 


24Mg 


0.39 


0° 


0.0 0.0 


14.12 9.51 


26Mg 


0.20 54° 


0.0 0.86 13.08 11.23 


28Mg 


0.0 


- 


0.0 1.03 


9.21 13.30 


30Mg 


0.0 


- 


1.31 1.03 


5.48 15.49 


32Mg 


0.0 


- 


0.0 1.03 


5.83 17.55 


34Mg 


0.37 


0° 


1.45 0.0 


4.12 20.18 


36Mg 


0.33 


0° 


1.43 0.0 


3.21 21.95 


38Mg 


0.30 


0° 


1.47 0.0 


2.38 23.69 


40Mg 


0.29 


0° 


0.91 0.0 


1.31 25.28 



but with the density- dependent contact interaction for the pairing channel. The space is 
truncated by the quasi-particle energy cut-off of Ecut = 60 MeV and also by the cut-off for 
the magnetic quantum number of the quasi-particle angular momentum, Qc = 19/2. In this 
QRPA calculation, the residual spin-orbit and Coulomb interactions are neglected. Thus, in 
order to make a comparison more meaningful, we also neglect the time-dependence of these 
potentials in the Hamiltonian h{t), during the time evolution. 

In panels (a) and (c) of Fig. [H we show the isoscalar quadrupole strength distributions 
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5 10 15 20 25 30 35 5 10 15 20 25 30 35 40 

E [ MeV ] 

FIG. 1: Calculated isoscalar quadrupole strength distribution for '^^Mg: (a) Cb-TDHFB with time- 
independent spin-orbit and Coulomb potentials, (b) Cb-TDHFB, (c) Deformed QRPA without the 



23]. The 



residual spin-orbit and Coulomb interactions [22]. and (d) Deformed QRPA calculation 
smoothing parameter of F = 1 MeV is used. 

{K = and 2) for ^^Mg, calculated with (a) Cb-TDHFB and (c) QRPA. The ground state 
has an axially symmetric prolate shape with finite pairing gaps for neutrons (Table [T]). The 
HFB calculation with the contact pairing interaction for the panel (c) produces a deformation 
and average neutron pairing gap of /3 = 0.37 and A„ = 1.7 MeV, respectively. Note that a 
renormalization factor, which was used in Ref. 22|, is set to be unity in the present QRPA 



calculation. The peak energies in these calculations are approximately identical, however, 
the height of the lowest peak is noticeably different. We suppose that this is due to the 
difference in the pairing energy functionals. 

In panels (b) and (d) of Fig. [H we show another comparison between the Cb-TDHFB 
calculation and the QRPA calculation of Losa et al[23] using the transformed harmonic 
oscillator basis. Since this QRPA calculation includes all the residual interactions, it is 
compared with the Cb-TDHFB calculation with the fully self-consistent time dependence. 
It turns out that the residual spin-orbit and Coulomb interactions slightly shift the giant 
quadrupole resonance higher in energy, while they shift the lowest peak lower in energy. 
Actually, these shifts are mainly attributed to the residual spin-orbit interaction and the 
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effect of the residual Coulomb is very small. The results in panels (b) and (d) well agree with 
each other, except for the height of the lowest peak. Again, this may be due to the different 
treatment of the pairing, because Ref. [23i] also uses the contact pairing interaction. These 
comparisons indicate that the small-amplitude Cb-TDHFB calculation well reproduces a 
fully self-consistent QRPA calculations. We would like to mention that, for the isovector 
dipole excitations, the agreement is even better than the isoscalar quadrupole cases. 



B. Isovector (El) dipole excitations 



Cb-TDHFB 

FAM 







FIG. 2: El strength distribution for 



Here, we discuss properties of the isovector dipole 
excitations including low-energy pygmy dipole res- 
onances (PDR) and high-energy giant dipole reso- 
nances (GDR). First, let us show in Fig. |2] that 
the comparison between results of the present Cb- 
TDHFB calculation and those of the RPA calcula- 
tion. The fully self-consistent RPA calculation has 
been done with the finite amplitude method (FAM) 



24Mg calculated with the Cb-TDHFB 
(solid line) and with the FAM[]J, UM 
(symbols). The smoothing parameter 
of r = 1 MeV is used. 
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24j . The same Skyrme func- 



developed in Refs. 

tional (SkM*) and the same model space have been 
used in these calculations. Since the ground state of 
the ^^Mg nucleus is in the normal phase (A = 0), 
these two results should be identical. This can be 

confirmed in Fig. |2l which clearly demonstrates the accuracy of our real-time method. 
Next, in Fig. |3l we show calculated El strength distribution for Ne isotopes. Here, 

S{E; El) is defined as 



S{E;E1)= Y. S,{E-E1)= Y, 5^|(n|F,|0)p(5(E-E„) 



(76) 



i=x,y,z 



t=x,y,z n 



where the one-body operator Fi is given by Eq. (!67|) . The K = strength is Sz{E;El) 
and K = 1 corresponds to Sx{E;El) + Sy{E]El). Here, for axially symmetric nuclei, the 
symmetry axis is chosen as 2;-axis. The ground states of ^"^'^^Ne are deformed in a prolate 
shape with A = for both protons and neutrons. Thus, the calculation is identical to the 
small-amplitude TDHF. Both nuclei show a prominent double peak structure. The lower 
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peak is located around 16 MeV and the higher one around 22 MeV. This conies from the 
deformation sphtting, and the lower peak is characterized as i^ = and the higher one as 
K = 1. The similar structure is seen in the neutron-rich nucleus ^^Ne. However, despite of 
the fact that the magnitude of deformation is roughly same as that of ^°'^^Ne, the position 
of the higher peak {K = 1) is lowered and the splitting is not as prominent as that in 
2o,22]\jg_ jj-^ oblate nuclei such as ^^Ne, the deformation splitting is not clearly seen in the 
total strength distribution, S{E]E1), because the high-energy peak becomes much smaller 
than the lower peak. 

For ^^~^^Ne, calculated ground states are in the superfluid phase for either neutrons or 
protons or both. Peak energies of the giant dipole resonance (GDR) gradually decrease as 
the neutron number increases, from about 20 MeV to 17 MeV. We have confirmed that 
the pairing correlation does not significantly affect the El strength distribution. However, 
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FIG. 3: Calculated El strength distribution for Ne isotopes. For deformed nuclei, the total strength 
([76D is decomposed into S;,{E; {El)) (thin solid line) and S^{E; {El)) + Sy{E; {El)) (dashed line). 
The z-axis is the symmetry axis for axially deformed cases. The smoothing parameter of F = 1 
MeV is used. 
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FIG. 4: Same as Fig. [3] but for Mg isotopes. 

for some cases, the ground-state deformation is changed by the presence of the pairing. 
For instance, the ^^Ne nucleus is deformed in the prolate shape if we neglect the pairing 
correlations. In contrast, the present BCS calculation produces the spherical ground state. 
The low-energy El strength, which is often called "pygmy resonance", is of significant 
interests. In Ne isotopes, there are two effects to create the low-energy El strength: One is 
a large deformation splitting which brings the lower peak down to around 15 MeV. Another 
effect comes from the neutron excess. In ^^~^^Ne, the pygmy peaks appear below 10 MeV. 
For ^^Ne, this low-energy peak structure has been recently measured at RIKEN |25| . The 
calculated pygmy position is around 8 — 9 MeV, which agrees with experimental data^] and 



with the other QRPA calculations 
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26J . For nuclei with even more neutrons {A > 28) 



there appear a double-peak structure below 10 MeV. 

In Fig. m El strength distribution for Mg isotopes are displayed. ^^Mg is nearly oblate, 
but has a triaxial shape with 754°. The low energy peak at 18 MeV is prominent in this 
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nucleus. In ^^Mg and heavier isotopes, there are pygmy states below 10 MeV. As is in 
Ne isotopes, there appear a double-peak structure for A > 30, though the heights of these 
pygmy peaks in Mg are lower than those in Ne isotopes. 

In order to investigate how the low-energy 
pygmy strength changes as the neutron number 
increases, we define the low-energy El ratio by 
mi{Ec)/mi with E^ = 10 MeV, where 



w 



mi{E) 



E'S{E'; El)dE', (77) 



I 2 

S 



Mg 




J I I I I I L_ 
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and nil = ?tii(oo). This ratio is shown for cal- 
culated even-even Ne and Mg isotopes in Fig. |5l 
Both isotopes with A^ = 10 ~ 14 have very lit- FIG. 5: Ratio of low-energy El energy- 
tie El strength below 10 MeV. Then, the ratio, weighted sum value to the total sum value, 
mi{Ec)/mi, shows a rapid increase as a function as a function of neutron number. See text 
of neutron number. Especially, we can see abrupt for details. 

jumps between A^ = 14 and 16, and between A^ = 20 and 22. The first jump between A^ = 14 
and 16 seems to be due to occupation of neutron si/2 orbital. In contrast, the second jump 
between A^ = 20 and 22 may be due to the onset of the deformation and the neutron pairing. 
The low-energy strengths in Ne isotopes show roughly twice larger values compared with 
Mg nuclei with the same neutron numbers. This may be attributed to the difference in the 
separation energy (chemical potential). 

Finally, let us present photoabsorption 
cross sections in the GDR energy region 
(i^^ = 10 ~ 30 MeV), together with exper- 
imental data|27|. For ^^Mg, the peak en- 
ergies of the GDR are underestimated by 
about 3 MeV. This disagreement has been 




25 10 15 

E [MeV] 



FIG. 6: Photoabsorption cross sections for already found in Ref.^ for ^^Mg. The 
24,26Mg. Experimental data (symbols) are taken present calculation also indicates that this 



from Ref. 



The smoothing parameter of underestimation of the GDR peak energy is 



1 MeV is used for the calculations. 



also true for ^^Mg. The El strength distri- 
bution for ^^Mg is very similar to that in 
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Fig. 12 (bottom panel) in Ref. 23|. In light nuclei, the GDR energy is systematically un 



derestimated in most of the Skyrme functionals 2J], that seems to be true for nuclei with 
superfluidity. 

VII. CONCLUSION 

We have developed an approximate approach to the time- dependent Hartree-Fock- 
Bogoliubov (TDHFB) theory, using the canonical-basis representation for time evolution 
of the densities and pairing tensors. Although, in general, the pair potential is not in a diag- 
onal form in the canonical basis, if it is approximated in such a form, the TDHFB equations 
can be enormously simplified, to give a canonical-basis TDHFB (Cb-TDHFB) equations. In 
this paper, we have treated a full Skyrme functional for the particle-hole channel, however, 
used a simple schematic functional for the pairing channel. Since the schematic pairing 
functional violates the gauge invariance, it requires a special choice for the gauge condi- 
tion. The Cb-TDHFB equations contain the TDHF as a special case of absence of pairing 
correlations. Its static limit is identical to the HF-I-BCS approximation. We have shown 
that the equations possess many of desired properties analogous to the original TDHFB 
theory, including the average particle number and the average total energy. We have also 
investigated analytical properties of its small-amplitude limit and found that the Nambu- 
Goldstone modes correspond to zero-energy normal modes and are automatically separated 
from other finite-energy modes. 

We have developed a computational program for real-time propagation based on the Cb- 
TDHFB equations using the three-dimensional (3D) coordinate-space representation. To 
test the accuracy and validity of the present method, we have calculated the isoscalar 
quadrupole strength distribution in deformed ^^Mg with the small- amplitude real-time 
method and compared with deformed QRPA calculations with a standard diagonalization 



method 
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231]. Results well agree with each other, except for the quadrupole strength of 
the lowest state located around 3 MeV. This may be due to the difference of the pairing 
energy functional used in the Cb-TDHFB and QRPA calculations. 

Then, we have calculated the El strength distribution for even-even Ne and Mg isotopes 
systematically. The ground-state properties of these isotopes changes from one nucleus to 
another. For instance, there are a variety of shapes including spherical, prolate, oblate, and 
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triaxial deformations. The gap energies also significantly changes, depending on the particle 
number and deformation. The 3D representation allows us to treat all of these nuclei in 
a self-consistent and systematic manner. Typical deformation splitting of the giant dipole 
resonances (GDR) is predicted for prolate deformed nuclei, 20,22,32j\j-g ^^^ 24,34-40]y[g^ rjj^g 
neutron-rich deformed nuclei, such as ^^Ne and ^^~^°Mg, show a K = peak around 15 MeV 
and a significant strength in a low-energy tail at 5 ~ 10 MeV. The low-energy El pygmy 
strength is almost negligible for ^°~^^Ne and ^^~^^Mg, but suddenly starts increasing at the 
neutron number 16 and another jump at 22. This seems to be due to the occupation of the 
neutron si/2 orbital and the onset of the neutron pairing. The effect of the deformation also 
plays a role for the increase of the pygmy strength in low-energy region. These low-energy 
El strength is of significant interest in studies of the element-synthesis reactions in stars 
and in explosive environments. 

The Cb-TDHFB method is easily applicable to heavier systems. Its computational task is 
roughly same as that of the TDHF. Furthermore, adopting the pairing functional calculated 
from an interaction instead of the schematic one, it can be used for calculation of the heavy- 
ion collision dynamics beyond the TDHF, including dissipative dynamics induced by the 
pairing interaction. 

Appendix: Proof of properties of Cb-TDHFB with Eg 

Among the properties listed in Sec. lIVBt the conservation of the orthonormal property 
and that of the particle number are trivially identical to Sec. IIIII In the foUowings, we show 
a simple proof of the other properties. 



1. Average total energy conservation 

Using the Cb-TDHFB equations, it is easy to show that the time derivative of the 
schematic pairing functional of Eq. (!60ll . gives 

k,l>0 ^ ^ fc>0 

Only with the special choice of the gauge parameters, fl62|) . we observe the conservation of 
the total energy. 
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2. Stationary solution 

Following the arguments in Sec. IIIIC| it is easy to see that the stationary solution corre- 
sponds to the ordinary HF+BCS result, but only when we adopt the gauge fixing fl^^ . 

3. Small amplitude limit 

With use of the pairing functional ( 160|) . we can no longer assume the time-dependent 
phase factor of Eq. f H^ for Afc(t). Instead, we only extract the global phase related to the 
chemical potential from Kk{t) and Afc(t). 

\Mt)) = m + \SMt)), \Mt)) = 14) + \SMt)), (A.2) 

f^kit) = {4 + SKkit)}e-'''^\ Afc(t) = {A" + 5Afc(t)}e-2^^*, (A.3) 

Pkit) = pl + Spkit), h{t) = ho + 6h{t). (A.4) 

Using the gauge condition (162|) . we have the following equations for the small-amplitude 
limit of the Cb-TDHFB equations. 

^Yt\^Mt)) = (ho - el)\5Mt)) + (1 - \4)mSh{tMl), (k o k), (A.5) 

^Q^^Pkit) = Al6Kk{t) + nlJ^GkiSKiit) - c.c, (A.6) 

l>0 

+{2pl - l)Y,Gki5t^i{t) + 2/\l5pk{t). (A.7) 

Here, we used the equation, {4>l + 5(pk{t)\{hQ+5h{t))\(t)l + S<j)k{t)) = (^1 + {4>l\^K^)\4>l) , which 
can be verified because of the norm conservation. 

a. Translation and rotation 

The same argument as that of Sec. IIIID II leads to 

ia{ho - el)S\cg) + (1 - |0°)(0°|)za[5, /.o]|0^) = 0, (A.8) 



where we multiply the projection (1 — |0fc)(0°|) on both sides of Eq. (!52|) . Equation (JA.SP 
means that |50f ) = i(l — |0fc)(0fc|)*S'|0^) and 5hs = ^[5, /iq] correspond to a zero-energy 
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normal-mode solution for Eq. (lA.Sp . Note that \S(j)f) = iaS\(j)^) and \S(j)f ) produce the 
identical density fluctuation. 6pk = and Suk = also satisfy Eqs. ( 1A.6I) and ( ]A.7p . since 
{4>k\^hs\(pl) = i{4>l\[S, ho]\(f)l) = 0. Therefore, the translational and rotational modes appear 
as the zero-energy normal modes. 

b. Pairing rotation 

When the ground state is in the superfluid phase, the transformation e*^^ changes the 
phase of k^. 

5fi:f = e^''4 -4^ 2i9Kt (A.9) 

5p^ = 6hN = 0, M) = |(50f ) = 0. (A.IO) 

Using Eq. (HOl) . it is easy to see that these quantities correspond to a normal- mode solution 
with a; = for Eqs. (1A.5I) . (IA.6I) . and ( 1A.7I1 . Thus, the pairing rotational modes appear as 
the zero-energy modes. Again, without the gauge condition (!62|) . we would not have this 
property. 

c. Particle-particle (hole-hole) RPA 

In case of k° = 0, it is easy to see that the p-p (h-h) channels are decoupled from the p-h 
channels. The p-p and h-h dynamics are described by the following equations. 

t^5Kk{t) = 2{el-X)6Kk{t)±J2GkiSKi{t), (A.ll) 

l>0 

where the sign + (— ) is for hole (particle) orbitals. Again, introducing the forward and 
backward amplitudes in the same way as Sec. IIIID3t Eq. (jA.lip can be rewritten in a 
matrix form as 



2epOpp' — Gpp' Gph' 

Ghp' —2e^5hh' — Ghh' I \^h' I yO — 1 / \Zh 

where Zp = Xp [Zp = Yp) and Z^ = Y^ {Z^ = X^) for the p-p (h-h) channel. 





(A.12) 
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